Inelastic collapse in one-dimensional driven systems under gravity 
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We study the inelastic collapse in the one-dimensional iV-particle systems in the situation where 
the system is driven from below under the gravity. We investigate the hard-sphere limit of the 
inelastic soft-sphere systems by numerical simulations to find how the collision rate per particle 
n co n increases as a function of the elastic constant of the sphere k when the restitution coefficient e 
is kept constant. For the systems with large enough N > 20, we find three regimes in e depending 
on the behavior of n co \\ in the hard-sphere limit: (i) uncollapsing regime for 1 > e > e c i, where 
n co n converges to a finite value, (ii) logarithmically collapsing regime for e c i > e > e C 2, where n co ii 
diverges as n co il ~ log k, and (iii) power-law collapsing regime for e C 2 > e > 0, where n co ii diverges 
as n co n ~ k a with an exponent a that depends on N. The power-law collapsing regime shrinks 
as N decreases and seems not to exist for the system with N = 3 while, for large N, the size 
of the uncollapsing and the logarithmically collapsing regime decreases as e c \ ~ 1 — 2.6/N and 
e C 2 — 1 — 3.0/JV . We demonstrate that this difference between large and small systems exists 
already in the inelastic collapse without the external drive and the gravity. 

PACS numbers: 45.70.Mg, 45.50.-j 



I. INTRODUCTION 

The inelastic hard-sphere system is one of the simplest 
models of granular media. It consists of rigid spheres 
that interact with each other only through instantaneous 
inelastic collisions. Minimum ingredients of the granular 
systems are taken in this system, for which the efficient 
event-driven algorithms have been developed for molec- 
ular dynamics simulations [lHU as well as sophisticated 
kinetic theories for analytical study (see, for example, 
Ref. [|). 

With this idealization of the granular media, how- 
ever, it has been known that infinite number of collisions 
among a finite number of particles can occur in a finite 
length of time. This phenomenon is called inelastic col- 
lapse [5, 6]. Process of collisions involved in the inelastic 
collapse has been studied for one-dimensional (1-d) (Hill 
and two-dimensional (2-d) systems (7rjTl|. and the con- 
ditions for the inelastic collapse have been obtained in 
some situations 

One of the simplest cases is the freely cooling granular 
gas, in which the inelastic hard-sphere system develops 
freely without any external forces 0, El ■ In the 2-d sys- 
tems, it has been shown that the particles that partake 
in the inelastic collapse form a string like linear structure 
for the case of frictionless particles Q, while they form 
a string like zigzag pattern for the case of the frictional 
particles [l(| ■ Another simple case is a simple shear flow, 
where the collapsing particles have been shown to form 
a linear string structure typically oriented along the di- 
rection 45° from the flow direction fijl . 
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In the hard-sphere idealization, once the inelastic col- 
lapse occurs, the system cannot proceed further without 
additional assumptions for the dynamics, such as those 
in the contact dynamics [12, [l]| . A simple way to escape 
from this difficulty is to suppress the inelastic collapse 
by employing the velocity dependent restitution coeffi- 
cient that goes to one as the colliding velocity goes to 
zero [ill, [HI [26| . In actual systems with finite rigidity, 
the inelastic collapse should never occur. 

Although the inelastic collapse is a singular behavior in 
the idealized system of the infinitely hard spheres with 
a constant restitution coefficient, its relevance to some 
physical behaviors has been suggested. In flowing con- 
figurations, it has been demonstrated by numerical sim- 
ulations that there exists strong correlation between the 
force chain network and the chain like structure formed 
by particles that collide repeatedly with each other in 
the hopper flow [l?} . The inelastic collapse has been also 
discussed in connection with the formation of correlation 
in the shear flow [lj| . 

In order to study how the inelastic collapse affects sys- 
tem behaviors in physical situations, it is natural to in- 
vestigate the soft sphere system with finite rigidity and 
see how the inelastic collapse appears in the hard-sphere 
limit. If you take, however, a simple limit of the infinite 
elastic constant with finite dissipation parameters, the re- 
sulting restitution coefficient tends to one, therefore the 
inelastic collapse does not occur. Thus the pertinent hard 
sphere limit for this purpose is the limit of infinite clastic 
constant with keeping the resulting restitution coefficient 
constant by making the dissipation parameter infinite. 

Mitarai and Nakanishi studied such limit by examining 
the limiting behavior of the collision rate n co n for the 2-d 
gravitational flow [l9j . The hard-sphere limit was taken 
as the limit of the infinite elastic constant k with the 
restitution coefficient e being kept constant. They found 
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that n co u converges to a finite value in the collisional flow 
regime, while it diverges as n co \\ ~ fc Q as k — > oo in the 
frictional flow regime. The exponent a was estimated to 
be about 0.4 in their case, i.e. the 2-d gravitational flow 
on a flat slope with ten layers of particles and the resti- 
tution coefficient e = 0.7. More recently, Brewster et al. 
studied the three-dimensional gravitational flow and ob- 
tained a ~ 0.25 for the system with 90~100 layers of par- 
ticles on a rough slope and e = 0.88 '[Hi . Although the 
divergence of collision rate implies emergence of inelastic 
collapse in the hard-sphere limit, a simple consideration 
on exponentially decreasing collision time interval would 
give the logarithmic divergence, and the mechanism for 
the power law divergence has not been understood yet. 

Motivated by these findings of the power law diver- 
gence in the gravitational slope flow, in this paper we 
take a closer look at the problem in an even simpler sys- 
tem, namely, a 1-d inelastic particle system under the 
gravity with an external excitation from a bottom of the 
system. The external excitation at the bottom is sup- 
posed to mimic the excitation by random collisions of 
particles with the slope in the gravitational flow, and our 
1-d system is intended to capture only the particle motion 
perpendicular to the slope. By numerical simulations, we 
will show that even this simple model exhibits the power 
law divergence of the collision rate, n co n ~ k a . 

This paper is organized as follows. In Sec. II, we begin 
by introducing our model and describe a method to study 
the hard-sphere limit of soft spheres in our simulations. 
The system with N = 3 is analyzed to show the logarith- 
mic behavior n C oii ~ log k in the hard sphere limit. In 
Sec. Ill, after describing the simulation procedure, first 
we present the simulation results for the systems with 
small number of particles (JV = 3^6); Only the loga- 
rithmic behavior in the hard sphere limit is observed for 
N = 3 while the power law divergence regime appears for 
the larger system in the smaller e region. Then we show 
the simulation results of systems with large number of 
particles (N > 20) and demonstrate that there exist the 
three distinct regimes for the limiting behavior of n co \i in 
e. We discuss the origin of these limiting behaviors based 
on the simulation results of inelastic collapse in the 1-d 
free space. Summary and conclusion arc given in Sec. IV. 



II. ONE-DIMENSIONAL MODEL OF 
GRANULAR FLOW 

A. Model 

We consider the 1-d model by focusing the particle 
motion only perpendicular to the slope (see Fig.Q]). The 
particles are allowed to move only along the z-axis under 
the influence of gravity and the lowest particle is excited 
by the bottom floor. 

Let us consider N identical particles with mass m and 
diameter d. The particles are numbered from the bot- 
tom starting with i = 1, and can interact only with 



their adjacent particles through the soft-sphere interac- 
tion. The excitation by the random collision with the 
slope is represented by the thermal floor located at the 
bottom z — 0. When the lowest particle (i = 1) collides 
with the bottom, it comes off with a random velocity v 
by the Maxwcll-Boltzmann distribution 



p{v) 



exp 



2k B T Q 



(1) 



where Tq is temperature of the thermal floor and fee is 
the Boltzmann constant. 

The interaction between soft spheres is given by the 
so-called spring- dashpot model @, [H|. Let z% and Vi de- 
note the coordinate and the velocity of particle i, re- 
spectively, then the overlap between the adjoining two 
particles i and i + 1 is given by ir^i+i = d — (z^+i — Zj). 
The relative velocity between i and i + 1 is denoted as 

= Vi - v i+ i = dx iii+ i/dt. Then, the force 
exerted on particle i by particle i + 1 is given by 



/1 



-kx it i +1 - flTM i|i+ i (for > 0), 

(for x i>i+ i < 0). 



(2) 



The first term of Eq. © represents the elastic force by 
the Hookean law with the elastic constant k. The second 
term denotes the dissipative force proportional to the rel- 
ative velocity i>j,i+i, where D is the damping constant. 
The force acting on the particle i + 1 by the particle i is 
given by /i+i,, = —fi.i+i- Note that the dissipative force 
is discontinuous at a^j+i = 0. 

The equation of motion for the particle i is then given 

by 

m ~JT = ~ m 9 + fhi+i - fi-i,h (3) 
dt 

where g is the gravitational acceleration. For our linear 
force law of Eq. ©, the duration time of contact for a 
binary collision r c is constant and given by 

(4) 



y/{2k/m) - D 2 



B. Hard-sphere limit 

For the linear force law in Eq. ([5]), the restitution co- 
efficient e of a binary collision is given by 



Vj,i+l\ t=Tc 
Vi,i+l\ t=0 



exp (-Dt c ) 



(5) 



flow 




O 

o 
o 
o 



gravity 



|<— thermal floor 

2D model ID model 

FIG. 1: A schematic representation of relation between the 
granular slope flow and the 1-d driven system. 
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using the duration time r c of Eq. (|4]) . By solving this for 
D, we obtain 

ZiT (me) 2 
Y m ' ^ 2 + (lne) 2 ' W 

The hard-sphere limit is defined as the limit of k — >• oo 
with keeping e constant. Thus, in this limit D diverges 
as D oc fc 1 / 2 and r c goes to zero as r c cx k~ x / 2 . 



For iV = 4, the critical value can be evaluated as 
e c (4) = e™ all (2) because of the symmetry in the order 
of collisions. If such symmetry in the order of collision 
process is assumed for N > 4, one may obtain the rela- 
tion 

e c (2iV) = er U (N), (11) 

but numerical simulations have shown that the relation 
Eq. QTJ) is not valid for large N [H]. 



C. Inelastic collapse 

Bernu and Mazighi [j| have studied N inelastic hard 
spheres thrown against a wall in 1-d system and showed 
the inelastic collapse can occur if e is less than a criti- 
cal value e™ all (7V). They gave an analytical expression 
for e™ a11 (N) using the independent collision wave (ICW) 
model: 



wall 



(iV) = tan 2 



1 - 



N 



(7) 



This is exact for the case of N — 2 but is an approxi- 
mation for N > 2 because the model ignores interaction 
between collision waves. Using another model called the 
cushion model, McNamara and Young [f| have obtained 
an estimate for the minimum number of particles jV^ 3,11 
that is required for collapse when the restitution coeffi- 
cient is e: 



N, 



ail,-, = ln(4/(l - e)) 
1 — e 



(8) 



This result becomes exact in the limit e — > 1. Compari- 
son between the ICW model and the cushion model has 
been discussed in Refs. [6|,[22j, and the numerical simula- 
tions show that the former is more accurate for iV < 15 
while the latter is better for N > 15. In the large N 
limit, both of the models give e™ a11 — > 1, but the asymp- 
totic forms are 



wall 

-cJCW 



for the ICW model and 



1 



7T 



wall 
^c, cushion 



l-_ln(4A0(l 



(9) 



lnln4iV\ , . 

■w ' (10) 



for the cushion model. 

The inelastic collapse can also occur in free space in 
1-d system if e is less than a critical value e c (N). A 
schematic picture of the three-body inelastic collapse is 
given in Fig. [2j McNamara and Young have shown 
that the three-body inelastic collapse can occur if e < 
e c (3) = 7 — 4V3, by using the 3x3 matrix M that 
relates the final velocities v' = (vi,v' 2 ,v' 3 ) of the three 
particles after two successive collisions with their initial 
velocities v = (vx,v 2 , V3) as v' = Aiv. 



D. Asymptotic analysis for N = 3 

In this subsection, we examine the asymptotic behavior 
of the total number of collisions n to t in a single collapsing 
event in the hard-sphere limit for the cases of N = 3 and 
show that n to t ~ logfc in the k — > 00 limit. 

The system with A^ = 3 is the smallest one where 
the inelastic collapse can occur, since the floor provides 
a thermal drive, thus the inelastic collapse can happen 
only in sequence of collisions among particles. Then, we 
can argue the behavior of n to t by considering a collision 
process of the three-body inelastic collapse in the hard- 
sphere limit as shown in Fig. [5J In this case, the time 
between the (n — l)th collision and nth collision, tj!| , 
between the same pair of particles, say 1 and 2, behaves 
as 



(n) 

12 



(n-l) 



q n t 



(0) 
12 i 



(12) 



where q is a constant smaller than unity (See Appendix 
A). 

Now, let us discuss the case of the soft spheres with 
a finite k and e < e c (3). In this case, initial binary 
collisions can follow a sequence similar to the inelastic 
collapse, but eventually all of the three particles are in 
contact after a finite number of collisions, and then fly 
away from each other with very small relative velocities. 
We estimate the total number n to t of collisions before all 
three particles are in contact at the same time. Similar 
estimation has been done for the case of a single inelastic 
soft sphere bouncing on a floor to show n to t ~ log k as 
k —> 00 [l9| . The three-body collapse like collision pro- 
cess shows essentially the same behavior as we will show 
in the following. 

(n) 

First of all, the collision interval t 12 for the case of the 
soft-sphere system is given by 



t 



in) 
12 



(13) 




FIG. 2: Schematic diagram for particle motion in the three 
body inelastic collapse in the one dimensional free space. 
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with the correction term At 12 m comparison with 
Eq. (|T2"j) because the collision duration r c is finite. It 
can be shown (see Appendix B) that 



-/t c , 



(14) 



where / is positive and a function of e. Substituting 
Eq. |14|) into Eq. 1(15)1. we obtain 



^12^ ~ ?^12 ^ f T c 
n-1 

i=0 



y c i2 



1 



-fTc 



(15) 



The number of collisions ntot before all three particles 
are in contact at the same time is given by the smallest 
n that satisfies the condition < 2r c , because two 

successive collisions between particles 1 and 2 cannot be 
shorter than the twice of the duration time. Thus, by 
requiring the relation 



("tot) 

12 



2t> 



for ntot 3* 1 and substituting Eq. (|15j). we obtain 



1 E 12 



1 - q" 



1-9 



(16) 



(17) 



n tot diverges and q ntot goes to zero in the limit k — > 00 
because r c oc fc -1 / 2 and < q < 1. Thus, Eq. (fT7|) can 
be written as g"* *^ ~ afc -1 / 2 , where a is a coefficient 
that depends on e. Therefore, we obtain n to t 



n to t 



21ogg 



log k + const. 



(18) 



If we assume that the frequency r of such a three-body 
process is independent of k, the collision rate is n co \\ ~ 
rntot and thus n co n ~ logfc in the hard-sphere limit. 



III. SIMULATION RESULTS 

The main quantity studied in this paper is the colli- 
sion rate per particle n co \\ defined as the average number 
of collisions (including collisions with the floor) per par- 
ticle per unit time for various values of parameters, fc, 
e, TV and To. We carried out numerical simulations to 
investigate n co ii in the hard-sphere limit. After describ- 
ing simulation procedure, we present the results for small 
number of particles 3 < ./V < 6 first, and then for large 
number of particles TV > 20. We find qualitative differ- 
ence between the two cases. 



A. Simulation procedure 

Numerical simulations are performed using the second- 
order Rungc-Kutta method with the time step dt = 



Tc/100, where r c is the duration time of binary collision 
given by Eq.(j4|). All particles are initially placed in such a 
way that there is no overlap between particles and veloc- 
ities are given randomly. After waiting for a sufficiently 
long time for the system to go through an initial tran- 
sient, we start taking data for various quantities and their 
time average. 

For numerical data, wc employ the unit system where 
the particle mass m, the diameter d, and the gravitational 
acceleration g are unities, 

m = d = g = l. (19) 

We set temperature of the thermal floor fee To = 1 unless 
otherwise stated. For a given set of TV and e, we measure 
the collision rate n coU for k = 10 5 , 10 6 , 10 7 , 10 s , 10 9 , and 
10 10 . 

Each collision event between two particles or between 
a particle and the floor is defined by their contact. The 
collisions between particles last for some duration time 
and they are counted everytime colliding particles sepa- 
rate, while the collisions with the floor are assumed to be 
instantaneous. The total number of collisions N co ii in- 
cludes the collisions between particles and those between 
a particle and the floor, and the collision rate per particle 
n co n is defined by 



n coU - Noon/ (NT), 



(20) 



where T is the simulation time length. We set T = 10 4 
for TV < 50 and T = 10 3 for TV > 50. 



B. Small systems 

Let us first consider systems with small number of par- 
ticles, in which a series of collisions occurs in a simple 
manner. Figure [3] shows the collision rates per particles, 
n co ii as a function of k for various values of e on the sys- 
tem with N = 3 ~ 6. For the system of TV = 3 (Fig. |3(a)[ ), 
the logarithmic behavior of n co u is clearly observed for 
e < e c (3) = 7 — 4-^/3, as has been suggested from the 
analysis in Sec. IID. On the other hand, for e > e c (3), 
ri co n converges to a finite value as k becomes large. It 
should be noted, however, that n co ii increases faster than 
logfc for e = 0.0718- e c (3). 

For the systems with TV =4, 5, and 6, such a region 
where n co ii increases faster than log k extends toward the 
smaller e region than e c (TV), e < e c (TV), as is seen in 
Fig. |31 Here, e c (TV) represents the critical restitution 
coefficient of the inelastic collapse for the free TV particle 
system evaluated by the ICW model with Eq. ([7]) and 
Eq. CO]), 



e c (TV) = tan 2 



-fl-- 

4 V TV 



(21) 



which we expect accurate for small N. 

In Fig. SI n co ii is shown as a function of e for k — 10 5 — 
10 10 on the system with N = 3 — 6. One may notice 
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FIG. 3: Collision rate per particle n co il vs - log k for various 
values of the restitution coefficient e for the systems with 
N = 3 (a), 4 (b), 5 (c), and 6 (d). 

that the curves have irregular-looking fine structures. We 
confirmed, however, that their statistical errors are small 
enough and that these fine structures are reproducible if 
we change sequences of random numbers in the simula- 
tions. Some of larger structures coincide with the criti- 
cal restitution constants e c (n) with (n = 3, • • • , N + 1), 
which are shown by the vertical dotted lines in Fig. 2) In 
the case of N = 3 (Fig. 4(a) I, one can observe a sharp 
peak at e c (3), and the peak value of n co \\ increases faster 
than logfc as k increases. For e < e c (3), n co ii increases 
by a nearly constant when k becomes 10 times larger, 
which suggests that n co ii increases logarithmically as is 
discussed above. 

For N = 3 ~ 6, there arc a couple of features in 
common. First, a sharp peak appears at e c (n)(n = 
3, • • • , N — 1) and the peak at e c (3) is highest. Secondly, 
a dip appears at e slightly larger than e c (N). Our simu- 
lation results (not shown here) suggest that this dip still 
exists for N = 10, becomes unclear for N = 25, and 
completely disappears for N = 30. 




FIG. 4: Collision rate per particle n co \l vs. e, plotted for 
k = 10 5 ~ 10 10 for the systems with N = 3 (a), 4 (b), 5 
(c), and 6 (d). The vertical dotted lines are at the restitution 
coefficient of e c (3), e c (4), ■ ■ • , e c (N + 1) from left to right. 

in the power-law for e < 0.8. The exponent a in the 
power-law regime depends on e, but is nearly constant 
a ~ 0.2 for e < 0.6 as can be seen from Fig. [He). In 
Fig. EJd), we also observed that the value of n co ii itself 
is nearly independent of e for any value of k for e < 0.6, 
where the exponent a is nearly constant. 

In the following, we will examine the transition region 
between the converging regime to the diverging regime 
carefully. Let us denote the lower limit of the restitution 
of the converging region by e c \ and the upper limit of 
the power-law diverging region by e C 2- A close look at 
the region 0.8 < e < 0.9 in the semi-logarithmic plots 
of Fig. EJa, b) reveals that there are two regimes within 
the region where n co ii diverges: the convex regime and 
the concave regime as a function of log k. In the convex 
regime, n co ii diverges faster than log k, suggesting that it 
is a part of the power-law regime. In the concave regime, 
the divergence is slower and it seems that n co \\ eventually 
shows the logarithmic divergence 



C. Large systems 

1. Collision rate in the large k limit 

For large systems, the power-law divergence of the col- 
lision rate dominates, but we can see clearly that there 
exists the region of restitution coefficient where n co n di- 
verges definitely slower than the power law. 

In Fig. O we plot n co ii for TV = 25 as a function of k for 
various values of e (a, b, and c), and as a function of e for 
various values of k (d). It is clear in the logarithmic plot 
of Fig. [SJc) that n co ii converges for e > 0.9 and diverges 



"■coil ~ b(e, N) log k + const. (22) 

in the large k limit with the coefficient b that depends on 
e and also on N. The lower limit of the converging regime 
e c \ is determined as the upper limit of the divergence; the 
data are fitted to Eq . ([22]) in the asymptotic region, then 
e c i is the point where 6 = 0. On the other hand, the 
value of e C 2 is estimated by the boundary between the 
convex and the concave regime. By these procedures, we 
obtain that e c i ~ 0.894 and e C 2 is somewhere between 
0.878 and 0.884 for N = 25. The values of e cl and e c2 
are determined for several values of N, and plotted with 
error bars in the 1/N-(1 — e) plane in Fig. [BJ One can 
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FIG. 5: Collision rate per particle n co n for N = 25. n co il — 
log k is plotted for (a) e = 0.80 ~ 0.90 with an increment 0.01 
from top to bottom, (b) e = 0.87 ~ 0.90 with an increment 
0.002 from top to bottom. logn co ii — logk is plotted for (c) 
e = 0.1 ~ 1.0 with an increment 0.1. (d) n co ii — e is plotted for 
k = 10 5 ~ 10 from bottom to top with the inset that shows 
in close-up near the critical values e c i ~ 0.894 and e C 2 ~ 0.88. 




see that they fit very well to the lines 
(1 



2 A a 

TV ' V 



fit \ 



3.0 
~N 



■if- ( 23 ) 



Note that their functional form is the same with the 
asymptotic form of e™Ji w in Eq. © . 

We further examine b(e,N) in Eq. (|2"2"|) as a function 
of both e and N by simulation data. From Eq. (|2"3"|) . we 
expect that b{e, N) is expressed by a simple function of 
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FIG. 6: Bifurcation diagram for the three regimes. The 
critical values of e c i and e C 2 are plotted for N =20, 25, 30, 40, 
50, 100, and 150. The error bars show ambiguity of the results 
in the procedure. The dashed and the dashed dotted lines 
show the fitting lines for e c \ and e C 2 by (1 — e^l) = 2.6/N and 
(1 — e^a) = 3.0/iV, respectively. The shaded region represents 




the region where the partially condensed state appears; 
boundary is estimated at N =20, 25, 30, 50, 100 and 150. 



its 



(1 — e) — A/N with a constant A w 2.6. This is actually 
what we find in Fig. (Ha), where we plot b(e,N)/N 5 / 2 
against (1 — e) — 2.6/N for various values of N and e in 
the logarithmic scale. One can see that the data collapse 
on a straight line with the slope 2, which means 



b(e, N) - N 5 ' 2 ((1 - e) - 2.6/iv) 2 

from which we confirm the asymptotic form 

2.6 



e c i 



1 - 



N 



(24) 



(25) 



It should be noted that this result is very close to the 
critical values of e below which clustering starts in the 
1-d granular system driven by a vibrating bottom plate, 
i.e. e c w 1 - 2.5//Y [H and e c 1 - 2.6/N [H. 

Based upon above analysis, we conclude that the crit- 
ical values, e c i and e C 2, do not coincide; therefore, in ad- 
dition to the converging regime, there are two diverging 
regimes in e with respect to behavior of n co ii in the hard- 
sphere limit k — > oo: (i) uncollapsing regime, e > e c i, 
where n co ii converges to a constant value, (ii) logarith- 
mically collapsing regime, e c \ > e > e C 2, where n co ii 
diverges as n co ii ~ logfc, and (hi) power-law collapsing 
regime, e < e C 2, where n co ii diverges as n co ii ~ k a . 



2. Partially condensed state 

We observe a partial condensation near the bottom 
around a certain value of e. Figure E] shows the spatial 
variation of the positional fluctuation (a) and the kinetic 
energy (b) of each particle for various values of e for 
N = 2b; Fig. 8(a) shows the standard deviation <Tj of the 
position of the particle i divided by the one ooi fc> r the 
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FIG. 7: (a) b(e,N)/N 5/2 vs. (l-e)-2.6/JV in the logarithmic 
scale. b(e,N) is defined in Eq. (|22p and estimated from the 
plots similar to those in Fig[5|c) for various values of N and e 
in the logarithmically collapsing regime. The dashed line gives 
a fit by 1.5((1 - e) - 2.6/iV) 2 for b/N 5/2 . (b) b(e, N)/N 1/2 
vs. N(l — e) in the linear scale using the same data as in 
(a). The range of values of e used to plot (a) and (b) is the 
following: 0.835 < e < 0.900 for N = 20, 0.876 < e < 0.900 
for N = 25, 0.894 < e < 0.920 for N = 30, 0.918 < e < 0.940 
for N = 40, 0.932 < e < 0.960 for N = 50, 0.963 < e < 0.972 
for N = 100, and 0.972 < e < 0.979 for N = 150. 

clastic case e = 1, and (b) shows the kinetic energy Ki of 
the particle i. Note that Ki — (1/2)/cbTo for any particle 
when e = 1. 

For e > 0.91, Cj/(Xjo and -?Q are larger in the region 
closer to the bottom and they decrease monotonically as 
the particle index i increases. This is because the ther- 
mal wall at the bottom supplies the kinetic energy to the 
bottom particle, and the kinetic energy is dissipated as 
it is transported away from the bottom via the inelas- 
tic collisions. However, around e ~ 0.9, a dip appears 
near the bottom both in o^/aoi and Ki, and there ap- 
pears the inversion layer where the temperature increases 
with i. This means that the low temperature and high 
density domain appears near the bottom. We call this 
the partially condensed state. For N = 25, the value of 
e ~ 0.9 where the partially condensed state appears al- 
most coincides with, but seems to be slightly larger than 
e c i — 0.894, i.e. the critical point of the inelastic collapse. 

The condensed domain with low Oi/oiQ in Fig. 8(a) 
extends towards the upper part of the system as e is 
decreased down to 0.84, where the whole system is con- 
densed. That is, the partially condensed state appears 
for 0.84 < e < 0.90, namely, the partially condensed 
state appears both in the logarithmically collapsing and 
the power-law collapsing regimes. 

In Figure [51 the region for the partially condensed state 
is shown by the shaded area in the 1/N-(1— e) plane. One 
can see that the upper bound of e (the lower boundary in 
Fig. IH]) for the partially condensed state nearly coincides 
with e c i for all of the cases examined. 
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FIG. 8: Spatial variation of the particle fluctuation of posi- 
tion and the kinetic energy for e = 0.82 ~ 0.92 for the system 
with N = 25. (a) The standard deviation of particle posi- 
tion <Tj of particle i, normalized by the corresponding value 
for e — 1, ao,i- (b) The kinetic energy Ki of the particle i. 
The data are shown for the cases with the elastic constant 
k = 10 10 . 

all the cases. In Fig. [9J the exponents a for e =0.1, 0.2, 
0.3 and 0.4 are plotted against logiV. One can see that a 
does not change by e but depends linearly on log N and 
can be fitted to 

a Rt = 0. 18 log N- 0.054 ~ 0.18 log (N/2). (26) 

This logarithmic dependence of the exponent a on N 
means that n co ii is given by 



"coll 



(N/2) 



0.181 



(27) 



for e < 0.4. 



4- Effect of the floor temperature Tq 



3. Exponent in the power-law collapsing regime 

We observed that n co ii is almost constant when e < 0.4 
for N — 25 as can be seen in Fig [5] We examined this for 
the case with TV = 20 ~ 150, and found this is true for 



We find no qualitative difference in the fc-dependence 
of "coil by changing ksTo from 1 to 10 in both systems 
with N — 25 and 50. We show in Table U the critical 
values e c i and e c i , and the average exponent of the power 
law a for fc B T = 1 (see Sec. IIIC.l), 2, 4, and 10 for the 
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FIG. 9: Exponent of the power law a for e = 0.1, 0.2, 0.3, 
and 0.4 plotted as a function of TV. The dashed line gives a 
fit by 0.18 log N - 0.054. 



systems with TV = 25 and 50. Here a is the arithmetic 
mean of the four values of a at e = 0.1, 0.2, 0.3, and 0.4, 
for which a's are almost constant independent of e. The 
both critical values and the power law exponent seem to 
be independent of T . 



TV 


A:bTo 


e c i 


e C 2 




a 


25 


1 


0.894 


0.878 ~ 


0.884 


0.20 




2 


0.892 


0.874 ~ 


0.880 


0.20 




4 


0.892 


0.872 ~ 


0.880 


0.20 




10 


0.890 


0.872 ~ 


0.878 


0.19 


50 


1 


0.944 


0.942 ~ 


0.936 


0.25 




2 


0.944 


0.938 ~ 


0.942 


0.25 




4 


0.944 


0.938 ~ 


0.940 


0.24 




10 


0.944 


0.934 ~ 


0.944 


0.24 



TABLE I: The critical values e c i, e C 2 and the average expo- 
nent of the power law a for various values of the floor temper- 
ature To for the systems with TV = 25 and 50. The average 
exponent a is defined as the arithmetic mean of the four val- 
ues of a at e = 0.1, 0.2, 0.3, and 0.4. 



D. Inelastic collapse in the free space 

In order to narrow down the possible origin of the 
power law divergence of the collision rate, we further 
simplify the system and consider the TV-particle system 
in the 1-d free space without the external drive and the 
gravity. We performed MD simulations to see how the to- 
tal number of collisions behaves in the hard sphere limit. 

In the initial state, TV particles of the diameter d are 
placed at an equal interval with the space a, 



Xi = \ i - 



TV + 1 



1,2, 



■N, 



with the initial velocities, 

Vi = -v sgn(x i ) + 8v &, 



(28) 



(29) 



where sgn(a;) is the sign function and £j's are random 
numbers distributed uniformly over the interval [—1,1); 
i>o and 5vq are positive parameters. 
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FIG. 10: ntot vs. log k for various values of e (a) for TV = 
3 in the log-linear scale and (b) for TV = 25 in the log-log 
scale. The parameters for the initial state are a = vq = 1 and 
5vo = 0.1. Each data point represents an average over 1000 
realizations by different random number £j's. The dashed 
lines in (a) show the asymptotic behavior given by Eq.{TH]) 
with q given by Eq. ()A12|) for the corresponding e values after 
the constants being adjusted to the data. The dashed line in 
(b) is the line with the slope 0.2. 



We count the total number of collisions n to t until the 
relative velocity of the end particles Vn — V\ becomes 
positive. The results are shown in Fig[TU] for the systems 
with TV = 3 (a) and 25 (b) for d = a = vq = 1 and 5vq = 
0.1. One can see that the total number of collisions ntot 
behaves in an analogous way with the collision rate n c . 



in the driven system under the gravity shown in Figs 3(a) 
and[5jc). For the case of TV = 3, n to t converges to a finite 
value when e > e c (3) ~ 0.0718 and diverges as log k when 



e < e c (3). The dashed lines in Fig 10(a) show Eq. (fTg|) 
with q given by Eq. (|A12[) for the corresponding e values, 
and adjusted constants. As for the case of TV = 25, n to t 
diverges as k a when e < 0.6 with the exponent close to 
the value a ~ 0.2 for the previous case with the drive 
and the gravity. 
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IV. SUMMARY AND DISCUSSIONS 

We have studied the inelastic collapse in the 1-d sys- 
tem under the gravity with the random driving from 
the bottom floor. Using MD simulations for the soft- 
sphere systems, we calculated the collision rate per par- 
ticle n co n and see how it diverges / converges in the hard- 
sphere limit; The hard-sphere limit is taken by the infi- 
nite limit of the elastic constant, k — > oo with the resti- 
tution coefficient e being kept constant. We have found 
that there are three regimes in the restitution coefficient 
e: (i) the uncollapsing regime for 1 > e > e c i, where 
ri co n converges, (ii) the logarithmically collapsing regime 
for e c i > e > e C 2, where n co n diverges as n co n ~ logfc, 
and (hi) the power-law collapsing regime for e C 2 > e > 0, 
where n co n ~ k a . For small TV systems, the region of e for 
the power-law collapsing regime is small and disappears 
for TV = 3. On the other hand, for large TV systems, the 
region for the power-law collapsing regime expands in the 
way that both of e c i and e C 2 approaches 1 as Eq. ((23")) . As 
for the floor temperature effect, we have checked the crit- 
ical restitution constants e c \ and e C 2 and the power law 
exponent a for the system of TV = 25 and 50, and found 
virtually no change for all of them in the temperature 
range of 1 < k B T < 10. 

If the intervals of collisions follow a geometrical se- 
quence toward the inelastic collapse, the logarithmic di- 
vergence of the collision number can be understood based 
on the consideration that the collision sequence termi- 
nates at the point where the collision interval becomes 
comparable with the duration time of binary collision. 
In the case of one particle bouncing on a floor under the 
gravity, it is obvious that the collision times follow the 
geometrical sequence. We have shown that it holds also 
for the three-particle system in the 1-d free space. Thus, 
the logarithmic divergence of the collision rate for the ex- 
ternally driven system can be understood if the inelastic 
collapses occur at a certain rate and they do not interfere 
each other nor are affected by the external drive in the 
hard-sphere limit, k — > oo. 

On the other hand, the power-law divergence of the 
collision rate is more intriguing. It has been reported in 
the gravitational slope flows [19|, but our results show 
that it occurs in an even simpler system, i.e. a 1-d ex- 
ternally driven system under the gravity. If we try to 
understand this in the same way as above, the collision 
interval should decrease as a power of collision number. 
This possibility is supported by the fact that the total 
number of collisions in the free space also diverges in the 
power law in the hard sphere limit. If this is true, we still 
need to understand how the power law sequence collision 
intervals arises. 

For the small TV systems, the collision rate shows a 
certain structure as a function of e at e = e c (n) of 3 < 
n < TV + 1, i.e. the critical restitution coefficient for the 
n-particle system in the free space. For TV = 3, there is a 
sharp peak at e = e c (3) and a dip at e = e c (4), while for 
3 < TV < 6, we find peaks at e = e c (n) for 3 < n < TV — 1, 



a dip or a shoulder at e ~ e c (TV), and a somewhat broad 
peak around e ~ e c (N + 1). Such a structure becomes 
vague for larger TV. Note that the relative motion of 
n particles under gravity with respect to their center of 
mass is equivalent to the motion in the free space as long 
as they do not interact with their surrounding particles. 
Thus, the structure at e c (n) for n < TV should be an effect 
of the inelastic collapse in which a part of the system is 
involved, but we do not understand yet how it shows up 
as a sharp peak or a dip structure, depending on the 
number of particles involved. 

For large TV, it is also intriguing that the collision rate 
is independent of e for a rather wide range; In the case 
of TV = 25, n co u is constant in the region e < 0.6 for any 
k, thus the exponent a also does not depend on e in the 
same region. Within the range 20 < TV < 150, 7i co n is 
almost constant for e < 0.4, but the power-law exponent 
increases depending linearly on log TV as TV is increased. 

One may find that some data points for TV > 100 in 
Figs. [6j \7\b), and [9] seem to deviate systematically from 
the asymptotic fitting expressions given by Eqs. (|23p . 

and (|2"6"]l . respectively. This may be due to the 
excessive load on particles near the floor. In the case of 
very large TV, the load on the particles near the floor be- 
comes so large that the contacts among them cannot be 
resolved into binary collisions but may remain as a long- 
lived contacts in the range of k investigated. In such a 
case, the system behavior may deviate from the assumed 
asymptotic forms. 

It is found that there appears the partially condensed 
state, where some particles near the bottom condense 
with lower kinetic energy. The region where the partially 
condensed state appears in the TV-(1 — e) plane covers the 
logarithmically collapsing regime; It starts hardly inside 
the uncollapsing regime and extends somewhat into the 
power-law collapsing regime. The condensed state near 
e c i contains only a few particles, but the number of con- 
densed particles increases as e decreases. The fact that 
the partially condensed state starts almost at e c \ suggests 
that the inelastic collapse causes the condensation, but 
it remains to be understood how a small number of par- 
ticles can condense by the inelastic collapse at e w e c i, 
that is much larger than e c (n) for small n. 

The partially condensed state has already been ob- 
served in the 1-d granular systems driven by a vibrating 
bottom plate in various forms of vibration. In the case 
of the sinusoidal vibration, the condensed state has been 
shown to appear for TV(1 — e) > 2.5 [23l] . and in the cases 
of a sawtooth vibration and a piecewise quadratic vibra- 
tion, for TV(1 - e) > 2.6 (24] • Our result shows that the 
partially condensed state appears below e c i, which means 
— e ) ~ 2.6 from Eq. ([2~5j) . These results are consis- 
tent with each other and show that the point where the 
system starts to condense is not sensitive to the driving 
mode. 

In summary, we have demonstrated that the inelas- 
tic collapse shows up in the 1-d driven system under 
the gravity as the diverging collision rate in the large k 
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limit with keeping the restitution coefficient e constant. 
By numerical simulations, we found that there are three 
regimes for the way that the collision rate diverges, i.e. 
the uncollapsing regime, the logarithmically collapsing 
regime, and the power-law collapsing regime. 



This relation can be expressed using the ratio of the rel- 
ative velocities to'™' = U32 /^l and m' = f'32 /f'21 
as 



(n) 



1 m^W™" 1 ) („_!) 
i* (n/™- 1 )) 2 ^ ' 



(A5) 



Appendix A 

In this appendix, we consider the three-body inelastic 
collapse of the inelastic hard spheres in the free space and 
derive the asymptotic behavior Eq. (fT2f of the time 
between two successive collisions between the particles 1 
and 2. 

Let v- n ' be the time of the nth collision between the 
particles 2 and 3 and t'^ be the time of the nth collision 



(n) 



r 



(n) 



between the particles 1 and 2, and define t\ 
iOO and = i (n+1) - t' {n) (See, Fig. HTJ). Similarly, 
the particle velocities just after t( n ) are (i = 1, 2, 3), 
those just after t'^ n ' are v'^ (i = 1,2,3). The relative 



velocities are denoted by = v 2 ' 



,(n) 



(n) 



1 and u 32 



(n) 



The separation between the particles 1 and 2 
at i'™) is denoted by and that between the particles 
2 and 3 at t' (n) by ar' 
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FIG. 11: Schematic picture of collision sequence occurring in 
the three-body inelastic collapse in the free space. 



The time t{ n ^ and to' 1 ' are then written as 



,(n) 



An) 



(n) 
C 21 



„/(«) 
6 32 

42 1 



respectively. The separations x^i and can 
expressed as 



(Al) 
also be 



(n) _ /(n-l) .(n-l) /(n) _ (n) ,(n) 

x 21 — y 21 ''2 1 32 ~~ y 32 n > 1-"-^) 

we can 



(A3) 



respectively. Combining Eqs. (|A1[) and ()A2j) . 
write 

„,/(«-!) 
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further rewrite it as 
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The sequence of to(") and m 1 ^, which are completely 
determined by the collision laws and the initial condi- 
tion m/ ) (or to'*'™'), has been studied by Constantin ct 
al. 25|. We briefly summarize their results that are rel- 
evant for our purpose in this appendix. Upon a collision 
between particles 1 and 2, velocities after the collision 
(v[, v 2 ) and before the collision (vi, v 2 ) are related by 



l-e l+e 
lie ' 



l 2 -e 



(A6) 



From this collision law, we can deduce the following re- 
lations: 



,(n) 



/("-!) 



/(«) 



bm 



/(n-l) ' 
b), 



(A7) 
(A8) 

where b = (1 + e)/2. If e < e c (3) = 7 - 4V3 and m (0) 
(or m ) is such that the collision sequence continues 
infinitely, i.e. the inelastic collapse occurs, then my 1 ' and 
m'W should converge to the stable fixed point values 



(-6 + Vb 2 - 4e) , 



2e 



-6- V^ 2 - 4< 



(A9) 
(A10) 



which are real when e < e c (3). 

Therefore, if both n and n' are so large that rrv- n > 



m ( n ) ~ m * anc l m 

written as 
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(n-l) 



, Eq. (|A5j) can be 

2(n-n') 
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It is straightforward to show 
2 ^ - ■ 2 



1 - 6e + e 2 - (1 + e)-/l - lie" 



1 - 6e + e 2 + (1 + e)Vl - 14e + e 2 



(A12) 

and < q < 1 for < e < e c (3). If we start to count the 
number of collisions at n', namely we put n' = 0, then 
we have 



t 



(n) 



(0) 



(A13) 



Because of symmetry with regard to exchange of par- 
ticles, t 2 should have the same property, and we finally 
obtain 



t 



(n) _ ^(n) _j_ ^(n-1) 



12 



<Z"t (0) 



12 J 



(A14) 



which is Eq. (IT2l . 
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t S(n) 



t S(n) 
L 2 



is the relative velocity of the particle 3 with respect 
to the center of mass of the particles 1 and 2, 



(») 



(n) In) 



V, 



(») 



,(n) 
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(B4) 
(B5) 



The relation Eqs. (|A1|) and (|A2[) for the hard particles 
should be modified for the soft spheres as 



FIG. 12: Schematic picture of collision sequence occurring 
in the three-body inelastic collapse like collisions among the 
inelastic soft spheres in the free space. 

Appendix B 

In this appendix, we turn to the problem of the three- 
body inelastic collapse like collisions among the inelastic 
soft spheres and derive the asymptotic behavior Eqs. (| 13[) 

and (|T4|) of the time between the instants of the end 
of contact at two successive collisions between the soft 
particles 1 and 2. We put superscript S for quantities that 
arc defined for soft particles in this appendix, in order 
to distinguish them from the corresponding quantities 
defined for hard particles in Appendix A. 

Let £ s ( n ) and f s ( n ) be the times of the beginning and 
the end of contact at the nth collision between the par- 
ticles 2 and 3, respectively (See, Fig. [T"2"|) . Similarly, let 

i' S and t'^™' be the times of the beginning and the 
end of contact at the nth collision between the parti- 
cles 1 and 2. We define the time intervals during which 

the particles move freely as t^ n ^ = t' SI ' ^ — t s (") and 
t S(n) = £S(n+i) _ t /S(n)_ Using the duration time r c of 

contact for a binary collision (sec Eq. Q), can be 
represented as 



,S(n) _ ,S(n) , .S(n-l) , „ 
L 12 — L l ' b 2 ' Al ct 

because two collisions occur during the time 



(Bl) 



We denote the relative distance between the particles 1 
and 2 just before the nth collision between the particles 2 
and 3 as cell > an( ^ thst just after the collision as x\^\ 
Similarly, the relative distance between the particles 2 
and 3 just before the nth collision between the particles 



1 and 2 is 



-S(n) 



32 



and that just after the collision is x 



/S(n) 
32 



The relative distances just before and just after a collision 
are related as follows: 
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x 32 — u 32 L l i 
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respectively. Combining Eqs. (|B2|) . (|B3|). (jB6j) and (|B7I) . 
we can write 
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Substituting Eqs. (|B4[) and ([B5]) into Eq. (JB8J) and using 
the ratio of the relative velocities and to'*-™-* defined 
in Appendix A, + r c can be expressed as 



,S(n) 



1 m^W™- 1 ) f t S(n-l) 
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(em 
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Note that the sequence of the particle velocities v^' (i = 

1, 2, 3) and that of and m'^ n ' are completely deter- 
mined by the collision laws and their initial conditions re- 
gardless whether the particles are hard or soft. Using the 
relations Eqs. flX7|) and (JMJ and the fact m'"', m' (n) < 
in the collapse like collision processes, it can be shown 
that the second term on the right-hand side of Eq. (|B9[) 
is negative. 

If m^™) and m converge sufficiently fast to their sta- 
ble fixed point values m* and m'*, and we start to count 
the number of collisions after w m* and m! ~ m'* 
are reached, we can write 



,S(n) 



(n) 

where V-{ is the relative velocity of the center of mass 
of the particles 2 and 3 with respect to the particle 1, and 
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C 21 
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32 
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m 
~~2 



{em 1 *) 2 



- 1 



,(B10) 



for any n > 1. 

Because of symmetry with regard to exchange of parti- 
cles, should have the same expression as Eq. (jBlOl) . 
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Substituting Eq. (|B10I) into Eq. (|Bl"j) . we finally obtain 
Eqs. |(I3) and USD: 



which is a positive function of e. 



,S(n) _ nf S(n-l) , 
''12 — 1 L 12 J 'cj 



(Bll) 



where 



/ = -m* 



(em'*) 2 



(B12) 
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